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Abstract. We propose a new iterative algorithm for generating a subset of 
eigenvalues and eigenvectors of large matrices which generalizes the method of 
optimal relaxations. We also give convergence criteria for the iterative process, 
investigate its efficiency by evaluating computer storage and time requirements 
and by few numerical tests. 
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The increasing computational power has stimulated a growing interest toward 
developing and refining methods which allow to determine selected eigenvalues of 
a complex quantum system with extreme accuracy. Widely adopted, specially for 
computing ground state properties, are the quantum Monte Carlo methods 0|, where 
a properly defined function of the Hamiltonian is used as a stochastic matrix which 
guides a Markov process to sample the basis. 

Alternatively, one may resort to direct diagonalization methods, like the Lanczos 
and Davidson |^ algorithms, of wide use in several branches of physics. The critical 
points of direct diagonalization methods are the amount of memory needed and the 
time spent in the diagonalization process. Because of these limitations, several systems 
are still out of reach even with the computer power now available. 

In this paper we present an iterative method, extremely easy to implement, for 
generating a subset of eigenvectors of a large matrix, give convergence criteria and 
show that it represents a generalization of the method of optimal relaxations . 

We assume first that the matrix A represents a self-adjoint operator A in an 
orthonormal basis {| 1), | 2), . . . , | N)} and is symmetric (a^- = {i \ A \ j) = aji). For 
the sake of simplicity, we illustrate the procedure for a one-dimensional eigenspace. 
The algorithm consists of a first approximation loop and subsequent iterations of 
refinement loops. The first loop goes through the following steps: 

la) Start with the first two basis vectors and diagonalize the matrix 



ai2 022 



where X\ = an. 



lb) Select the eigenvalue A2"^' and the corresponding eigenvector | 02^^) = K2^l 



for j = 3, . . . , iV 



Ic) compute foj^-* 



Id) Diagonalize the matrix 




(1) 



le) Select the eigenvalue A^^^ and the corresponding eigenvector | (ffj^^)- 
end j 

The first loop yields an approximate eigenvalue a|^^ = E^^"^ = X^q^ and an approximate 

eigenvector | -0'^^) =1 '/'tv^) =1 '/'o^') = J2iLi I *)■ With these new entries we start 
an iterative procedure which goes through the following refinement loops: 
for n = 2, 3, . . . till convergence 
for j = l,2,...,iV 

2a) Compute fej"^ = (</)5"\ \A\j). 

2b) Solve the eigenvalue problem in the general form 



det 




'33 




where the appearance of the metric matrix follows from the non orthogonality of the 
redifined basis | and | j). 

2c) Select the eigenvalue Aj""* and the corresponding eigenvector | 0^"''). 
end j 
end n. 
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The n-th loop yields an approximate eigenvalue 
eigenvector, at any step of the j-loop, we have 



a[,"+^\ As for the 



(n) 



q] I j) 



(1) 



with the appropriate normalization condition [p^"^]'^ 4- 
1. The iteration of Eq. (|^) yields the n-th eigenvector 



(n)l2 



2 n^"^^^"^/^^"^ 



AT 

1=1 



(2) 



where the numbers P- 



(n) 



p(n) _ 



are 

N 



(z = 0,l,...,7V-l) 



N 



1. 



(3) 



n 

The algorithm defines therefore the sequence of vectors (2), whose convergence 
properties we can now examine. The gj"^ and p^^^ coefficients can be expressed as 



(«) 
9 = 



B 



in) 



(«) 

J („): 



2 + 2if i-nK^i 



(4) 



where 



B 



(n) 



i-1 



A 



(») 



(") 



It is apparent from these relations that, if 

0, Vj-, 



(5) 
(6) 



the sequence | -0^"^) has a limit | which is an eigenvector of the matrix A. In fact, 
defining the residual vectors 

a direct computation gives for their components 



(7) 



^(") _ „(») 



{an - aJ"))(aS - A("))] " + gj;') {a,^. - A^")<5,a.} 



(8) 



In virtue of (g), the norm of the n-th residual vector converges to zero, namely 
II r^"^ ||— > 0. Equation (|^) gives therefore a necessary condition for the convergence 
of I '0'"'') to an eigenvector | tp) of A, with a corresponding eigenvalue E = lim i?'"'. 
This condition holds independently of the prescription adopted for selecting the 
eigensolution. Indeed, we never had to specify the selection rule in steps lb), le) 
and 2c). Equation (^) is not only a necessary but also a sufficient condition for the 
convergence to the lowest or the highest eigenvalue of A. In fact, the sequence A^"^ is 
monotonic (decreasing or increasing, respectively), bounded from below or from above 
by the trace and therefore convergent. 
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The just outlined algorithm has a variational foundation. Its variational 



counterpart is just the method of optimal relaxation M. Indeed, for the p^-"' and 



gj"'' given by Eqs. the a^"''(= (7j"V?'j"'') derivative of the Rayleigh quotient 
vanishes identically. 

On the other hand the present matrix formulation allows in a straightforward way 
for the optimal relaxation of an arbitrary number t of coordinates, thereby improving 
the convergence rate of the procedure. We only need to turn the two-dimensional into 
a. (t+ l)-dimensional eigenvalue problem in steps Id) and 2b), compute t elements bj 
in steps Ic) and 2a), and accordingly redefine the j-loops. The current eigenvector is 

still defined by the iterative relation (a^"'' = Q^^^ /p'j^^) 



which automatically fulfils the extremal conditions 




k) , (10) 



^ p(4")j = 0, fc = j + l,...,j+i. (11) 



Moreover, the algorithm can be naturally extended to generate at once an arbitrary 
number m of lowest eigenstates. We have simply to replace the two-dimensional 
matrices with multidimensional ones having the following block structure: A m x m 
submatrix diagonal in the selected m eigenvalues, which replaces Aj"\, a m' x m' 
submatrix corresponding to ajj and two m x m' off-diagonal blocks replacing fe^"'' or 

(n) 

Kj_\ j . This new formulation amounts to an optimal relaxation method of several 
coordinates into a multidimensional subspace. It avoids therefore the use of deflation 
or shift techniques for the computation of higher eigenvalues and eigenvectors. 

It remains now to investigate the practical feasibility of the method. The main 
issues to be faced are the storage and time requirements. In the one-dimensional case, 
we need to store a single A^-dimensional vector (the eigenvector) . The time is mainly 
determined by the j loop. This requires N operations for implementing point 2a) plus 
/c ~ 15 remaining operations. Since n — 1, 2, ric and j = 1, 2, N, the algorithm 
requires altogether ndN'^ + kN) operations. It follows that, for large dimensional 
matrices, the number of operations grows like iV^. For sparse matrices with an average 
number L of non zero matrix elements, the required number of operations is nc{L+k)N 
and therefore grows linearly with N. In the multidimensional case we need to store 
m A^-dimensional vectors. If necessary, however, we can keep only one at a time and 
store the remaining m — 1 vectors in a secondary storage. This latter feature clearly 
shows that the algorithm lends itself to a straightforward parallelization. Also in the 
multidimensional case, the number of operations grows as UcmN'^. 

The algorithm has other remarkable properties: i) It works perfectly even in case 
of degeneracy of the eigenvalues, ii) The diagonalization of the submatrices of order 
m + m' insures the orthogonalization of the full N-dimensional eigenvectors at each 
step. Therefore, no ghost eigenvalues occur, iii) The range of validity of the algorithm 
can be easily enlarged if we remove some of the initial assumptions. Clearly, the 
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iterative procedure applies to a non-orthogonal basis. We simply need to substitute 
steps la) and Id) of the first loop with the appropriate generalized eigenvalue problem. 
It applies also to non symmetric matrices. We have only to update both right and left 
eigenvectors and perform steps Ic) and 2a) for both non-diagonal matrix elements. 

In order to test the efficiency and the convergence rate of the iterative procedure, 
we have applied the method to several examples. The first is a 5-point finite difference 
matrix arising from the two-dimensional Laplace equation This is a block- 

tridiagonal matrix of rib 6-dimensional blocks, whose eigenvalues are 

' V 2(nb + l) 2(6+1)7' ^ ' 

where i = l,2...,nfc and j = 1,2...,6. As in Q, we considered a block-matrix 
with nf, = 15 and h — 20. We have tested the one-dimensional as well as the 
multidimensional version of the algorithm. As shown in Figure the iterative 
procedure converges much faster in the multidimensional case. In fact, the convergence 
rate increases with the number v of generated eigenvalues and is considerably faster 
than in Lanczos. It is also to be stressed that our algorithm allows for an arbitrarily 
high accuracy, up to the machine precision limit. The method, specially in its 
multidimensional extension, is quite effective even if applied to the same matrix with 
rib = b so as to allow for degeneracy. For rib = b = 80, it yields the lowest seven roots, 
including two couples of degenerate eigenvalues, with an accuracy of 10~^^. 

A second example, still taken from ||^, is a one-dimensional biharmonic band 
matrix whose eigenvalues 

Afc ^ l&sin^ , , fc = l,...,iV (13) 
2(iV-f-l) ' ' ^ ' 

are small and poorly separated from each other. A similarly high density of levels 
occurs in the Anderson model of localization Because of this peculiarity, the 
limit of the machine precision is reached for a modest increase of the dimension iV 
of the matrix. Our method reproduces perfectly any number of eigenvalues with the 
required accuracy after a small number of iterations. In the specific example discussed 
in |5) (A^ = 20) we attained the highest accuracy after eight iterations, much less 
than all methods discussed there. We have checked that, unlike others, our method 
works without any modification even if we increase the dimension N up to the limit 
compatible with the machine precision. In this case the number of iterations needed 
increase by an order of magnitude, in any case, below 100. 

A third example is provided by a matrix with diagonal matrix elements an = 
— a and off-diagonal ones = —a or a^- = according that they fall within or 
outside a band of width 2L. Such a matrix simulates a pairing Hamiltonian relevant 
to many branches of physics. We have considered a matrix of dimension = 10* 
and half-band width L = 400. We found convergence after 28 iterations, reaching 
an accuracy of 10~* for the eigenvalues. The time required to compute the lowest 
eigenvalue through the one-dimensional algorithm is i = 18697 s for a workstation of 
500 MHz and 512 Mb of RAM. 

Finally, we generalize the latter example by considering a full matrix of 
dimension N = 10^ with matrix elements an — + (— l)''+-^a . Their 

alternating signs are also to be noticed, since they decrease somewhat the rate of 
convergence of the process. We reproduce the lowest eigenvalue with an accuracy of 
10-^ 10"^, lO"'^, 10"* after ric = 42, 70, 155, 330 iterations respectively. 
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In conclusion, the present diagonalization algorithm is a generalization of the 
variational optimal relaxation method and, on the ground of the examples discussed, 
appears to be more competitive than the methods currently adopted. It seems to be 
faster and to require a minimal amount of computer storage. It is extremely simple to 
be implemented and is robust, yielding always stable numerical solutions. Moreover, 
it is free of ghost eigenvalues. Because of these features, we are confident that it can 
be applied quite effectively to physical systems, like medium-light nuclei or quantum 
dots with few electrons. 
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Figure 1. Convergence rate of the present algorithm apphed to a finite difference 
matrix deduced from Laplace equation for different numbers u of generated 
eieenvalues. The data referring to the Lanczos convergenge rate are taken from 
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